close all;
clc,clear;
warning off;
%%  导入数据(时间序列数据一列)
data=xlsread("Multi-Feature_Regression_Data.csv");
%% 样本
nums = length(data);   %样本数
re=data;
%%  数据归一化 索引
X=re(:,1:end-1);
Y=re(:,end);
[x,psin]=mapminmax(X',0,1);   %归一化
[y,pcout]=mapminmax(Y',0,1);   %归一化

%% 划分训练集与测试集
num=size(re,1); %总样本数
% k=input('是否打乱样本(是 1，不是0):');
k=1;
if k==0
    state=1:num;
else
    state=randperm(num);%增加泛化
end
r=0.8; %训练集占比
trainnum=floor(num*r);%训练集总样本数
xtrain=x(:,state(1:trainnum));
ytrain=y(state(1:trainnum))';

xtext=x(:,state(trainnum+1:end));
ytext=y(state(trainnum+1:end))';

%% 适应格式
kinds=size(x,1);%特征数
xtrain=reshape(xtrain,kinds,1,1,trainnum);
xtext=reshape(xtext,kinds,1,1,num-trainnum);
%%  创建网格
layers=[imageInputLayer([kinds,1,1])            %输入层
    % 第一层卷积
    convolution2dLayer([3,1],8,'Padding','same')%卷积枝大小3*1,8个卷积核
    reluLayer                                   %Relu激活层
    maxPooling2dLayer([2,1],'Stride',[1,1])     %最大池化层 池化窗口[2,1]
     % 第二层卷积
    convolution2dLayer([3,1],8,'Padding','same')%卷积枝大小3*1,8个卷积核
    reluLayer                                   %Relu激活层
    maxPooling2dLayer([2,1],'Stride',[1,1])     %最大池化层 池化窗口[2,1]

    dropoutLayer(0.1)                           %Droput层
    fullyConnectedLayer(1)                      %全连接层
    regressionLayer                             %回归层
    ];

%% 设置训练参数
options=trainingOptions('adam' ,...                              %梯度下降算法
                        'MiniBatchSize',32,...                   %批大小,每次训练的样本个数
                        'MaxEpochs',100, ...                     %最大迭代次数
                        'InitialLearnRate',0.01, ...             %初始学习率
                        'LearnRateSchedule','piecewise', ...     %学习率下降
                        'LearnRateDropFactor',0.1,...            %学习率下降因子
                        'LearnRateDropPeriod',80, ...           %经过多少次训练后学习率为*0.1
                        'Shuffle','every-epoch', ...             %每次训练打乱数据集
                        'Plots','training-progress', ...         %画出曲线
                        'Verbose',true);                         %观察进度

%% 训练模型
net=trainNetwork(xtrain,ytrain,layers,options);

%% 仿真预测
re1=predict(net,xtrain);
re2=predict(net,xtext);

%% 反归一化（实际值）
Ytrain=Y(state(1:trainnum));
Ytext=Y(state(trainnum+1:end));

%预测值
pre1=mapminmax('reverse',re1,pcout);
pre2=mapminmax('reverse',re2,pcout);

%% 预测
newdata=readmatrix("new_data_for_prediction.xlsx");
newy=newpre(newdata,psin,pcout,net);
%作图
figure
plot(newy)
xlabel('样本点');
ylabel('预测值');
%% 均方根误差
error1= sqrt(sum((pre1-Ytrain).^2)./trainnum);
error2= sqrt(sum((pre2-Ytext).^2)./(num-trainnum));

%% 相关指标计算
%R^2
R1=1-norm(Ytrain-pre1)^2/norm(Ytrain-mean(Ytrain))^2;
R2=1-norm(Ytext-pre2)^2/norm(Ytext-mean(Ytext))^2;

% MAE
mae1=mean(abs(Ytrain-pre1));
mae2=mean(abs(Ytext-pre2));

%% 图
figure 
plot(1:trainnum,Ytrain,'r-^',1:trainnum,pre1,'b-^');
legend('真实值','预测值');
xlabel('样本点');
ylabel('预测值');
title('训练集预测结果对比');

figure 
plot(1:num-trainnum,Ytext,'r-^',1:num-trainnum,pre2,'b-^');
legend('真实值','预测值');
xlabel('样本点');
ylabel('预测值');
title('测试集预测结果对比');

%百分比误差图
figure 
plot((pre1-Ytrain)./Ytrain,'b-o','LineWidth',1);
legend('百分比误差');
xlabel('样本点');
ylabel('误差');
title('训练集百分比误差曲线');

figure 
plot((pre2-Ytext)./Ytext,'b-o','LineWidth',1);
legend('百分比误差');
xlabel('样本点');
ylabel('误差');
title('测试集百分比误差曲线');

%拟合图
figure
plotregression(Ytrain,pre1,'训练集', ...
    Ytext,pre2,'测试集');
set(gcf,'Toolbar','figure');
function y=newpre(newdata,psin,pcout,net)
     %归一化
     x=mapminmax('apply',newdata',psin);
     kinds=size(x,1);
     num=size(x,2);
     %适应格式
     xnew=reshape(x,kinds,1,1,num);
     %预测
     pre=predict(net,xnew);
     %反归一化
     y=mapminmax('reverse',pre,pcout);
end